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1. Introduction 


ABSTRACT 


The power management strategy (PMS) plays an important role in the optimum design and efficient 
utilization of hybrid energy systems. The power available from hybrid systems and the overall lifetime of 
system components are highly affected by PMS. This paper presents a novel method for the determina- 
tion of the optimum PMS of hybrid energy systems including various generators and storage units. The 
PMS optimization is integrated with the sizing procedure of the hybrid system. The method is tested ona 
system with several widely used generators in off-grid systems, including wind turbines, PV panels, fuel 
cells, electrolyzers, hydrogen tanks, batteries, and diesel generators. The aim of the optimization problem 
is to simultaneously minimize the overall cost of the system, unmet load, and fuel emission considering 
the uncertainties associated with renewable energy sources (RES). These uncertainties are modeled by 
using various possible scenarios for wind speed and solar irradiation based on Weibull and Beta proba- 
bility distribution functions (PDF), respectively. The differential evolution algorithm (DEA) accompanied 
with fuzzy technique is used to handle the mixed-integer nonlinear multi-objective optimization prob- 
lem. The optimum solution, including design parameters of system components and the monthly PMS 
parameters adapting climatic changes during a year, are obtained. Considering operating limitations of 
system devices, the parameters characterize the priority and share of each storage component for serving 
the deficit energy or storing surplus energy both resulted from the mismatch of power between load and 
generation. In order to have efficient power exploitation from RES, the optimum monthly tilt angles of 
PV panels and the optimum tower height for wind turbines are calculated. Numerical results are com- 
pared with the results of optimal sizing assuming pre-defined PMS without using the proposed power 
management optimization method. The comparative results present the efficacy and capability of the 
proposed method for hybrid energy systems. 

© 2011 Elsevier Ltd. All rights reserved. 


due to their direct dependence on climatic conditions. In addition, 
the variations of the wind and solar energy may not match with load 


Distributed generation (DG) has been recently nominated as 
one of the solutions to reliable, less costly, and more efficient 
energy supply systems. Specifically, small DG systems with power 
level ranging from 1kW to 10MW [1], located near the loads, 
are extensively utilized both in grid-connected and stand-alone 
modes. Among available DGs, renewable-based systems (RES) 
such as photovoltaics and wind turbines have attained the most 
popularity due to ever increasing concerns about depletion of 
fossil fuels and global warming. They have also been getting more 
cost-effective during recent years. 

However, a significant drawback associated with solar and wind 
energy systems is their intermittent and unpredictable behavior 
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changes [2]. The reliability of the system is important to both plan- 
ning and utilization stages. Designing energy systems including 
solar and wind energies together, to some extent, reduces the depth 
of the problem. Since, the advantage of one source can overcome 
the disadvantage of the other one and vice versa [3]. In addition, 
taking into account the intermittency and uncertainty associated 
with solar and wind energy sources, improves the adaptation of 
design results with practical and realistic conditions. Another con- 
ventional solution is appropriate incorporation of energy storage 
devices such as batteries and hydrogen storage systems compris- 
ing electrolyzers, hydrogen tanks and fuel cells into the system 
[3]. Hence, hybrid energy systems are becoming more attractive 
to power engineers. 

One of the most prominent issues regarding hybrid energy 
systems is to determine their optimum design and operation 
modes taking account of regional conditions and load demand 
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characteristics. Therefore, many efforts have been made for design- 
ing and planning hybrid DG systems [1-18]. These works can be 
reviewed from various points of view such as system modeling, the 
applied data in the solution algorithm, objective functions to be 
optimized and the mathematical tool used to handle the optimiza- 
tion problem. In [2], the overall cost of the introduced hybrid system 
throughout the estimated lifespan of the installation has been 
minimized. The reliability of serving energy and proper functioning 
of components has been involved in the applied algorithm as a con- 
stant constraint not to be violated. Solar radiation and wind speed 
are also assumed deterministic. Particle Swarm Optimization (PSO) 
algorithm has been applied to handle the single-objective problem. 
In [3], the cost analysis of a hybrid wind/PV/fuel cell system has 
been focused on through a residential consumer case study. 

In [4-10], the cost has been minimized by selecting the corre- 
sponding system configuration and/or operation strategy. The load 
demand has been supposed to be completely or partially fulfilled 
with a fixed user-defined value for unmet load. In this case, the 
solution method has no choice to reach a compromise between 
cost and unmet load. In [9,10], pollutant emissions have been rep- 
resented in the form of cost and limited by being economically 
involved in the cost function, and thus the results depend on the 
cost coefficient assigned to the emissions. More realistic optimal 
configuration may be achieved by involving pollutant emissions 
measurement in their relative unit (kg per liters of fuel) instead of 
converting it to the cost. From the aspect of the system modeling 
accuracy, Ref. [5], has taken a fixed PV panels tilt angle into account 
using Hay and Devis (1980) and Orgin and Holland (1977) meth- 
ods, whereas [6] has imported the angle value as an optimization 
variable in the solution algorithm. In the latter work, commercially 
available types for each device and their costs have been taken 
into account, and optimization results representing size and type 
of each device have been demonstrated. The load demand has been 
assumed to be entirely fulfilled. The maximum power point tracker 
(MPPT) has been included in both PV and wind turbines systems. 
The presented approach in [7] has considered the influence of the 
uncertainty in solar irradiation on the sizing process of an off-shore 
PV power system using three probabilistic models such as Markov 
Chain. A size optimization procedure has been demonstrated in 
[8], where the results have been compared by using the optimiza- 
tion software called HOMER [19]. It has been stated that because of 
the nonlinearity and complexity of hybrid systems, the application 
of evolutionary algorithms like genetic algorithm generates better 
results than the application of classical optimization methods [8]. 

On the case of PMS [5], as one of the earliest works, has 
introduced two parameters named SDM and SAR, representing 
the minimum and maximum State of Charge (SOC) of batter- 
ies, respectively. The optimum values of the SOCs have been 
calculated within the context of the size optimization problem. 
Although [6] has applied an accurate model for PV and wind 
turbines systems representing installation parameters, the oper- 
ation control has been fixedly defined by the user and did not 
conform with the system conditions. In [8], the effect of the min- 
imum and maximum limit of SOC on the system operation and 
cost has been evaluated. The algorithm suffers a computation bur- 
den since a sub-algorithm for finding optimum values of control 
parameters has been executed for any single control vector gener- 
ated in the main optimization algorithm as a nested optimization 
loop. 

In [9], a grid-connected hybrid energy system, capable of 
mutual exchange of power with the grid, comprising wind turbine, 
micro-turbine, diesel generator, photovoltaic array, fuel cell and 
battery storage has been analyzed and on-line optimal power 
management has been attained. The problem has been treated as 
a single objective problem considering all objectives such as fuel 
emissions in terms of cost. The fuel emission has been factorized 


into three main gas components and values for each gas has been 
separately presented. 

The cost minimization of a hybrid energy system including 
hydrogen storage has been presented in [10]. Using the classic 
linear programming optimization tool, a comparison between an 
optimized dispatch strategy and a fixed one has been performed so 
that the latter strategy can be reformed using the values obtained 
for the dispatch variables of the former one. Ref. [11] has incorpo- 
rated the effect of the number of battery charge/discharge cycles 
relative to their depth of discharge (DOD) on their lifespan, oper- 
ation strategy and consequently the system cost. A similar work 
has been addressed in [12], which has modeled the uncertainty 
involved in operating and design characteristics of the system tak- 
ing the advantage of Stochastic Annealing optimization algorithm. 
This work has referenced the classification of uncertainties pre- 
sented in [13] as model-inherent, system inherent and external 
uncertainties existing in the system. 

In [14], a method has been introduced, in which the optimized 
configuration and operation of the system has been achieved 
based on twelve parameters defined relative to operational limits 
of different storage devices. This work has been rearranged in 
[15] considering a more complex hybrid system and a couple of 
objective functions including cost, emission and unmet load. 

As seen in [14], the power available from such systems and 
the overall lifetime of system components are highly affected by 
the applied power management strategy; hence, various strate- 
gies result in various designs for the system aiming to meet load 
requirements within the estimated lifespan of the system. 

In this work, a general method for optimal PMS of autonomous 
systems including various generators and storage media, com- 
bined with optimal design of system components is proposed. 
To do so, appropriate model of system components for the plan- 
ning problem is reviewed. To solve the introduced problem, the 
differential evolution algorithm (DEA), which is nominated to be 
capable of solving various nonlinear optimization problems, is uti- 
lized. The algorithm is accompanied with the fuzzy technique to 
handle the multi-objective optimization problem in a time effi- 
cient manner. In this case, there is no need to find appropriate 
coefficients as penalty factors of constraints violations and as 
objective function weights. A novel method, involved in the size 
optimization problem, is proposed to obtain values for the param- 
eters of PMS by which the optimum performance and minimum 
cost, emission and unmet load are achieved. Considering operat- 
ing limitations of system devices, these parameters characterize 
the priority and share of each storage component for serving 
the deficit energy or storing surplus energy both resulted from 
the mismatch of power between load and RES. Optimal values 
for design parameters and PMS parameters are simultaneously 
attained. 

The RES uncertainties are applied to the optimization procedure 
by scenario generation based on Weibull [20] and Beta distribu- 
tion functions for wind speed and solar irradiation, respectively 
[21]. Numerical results, including type and number of each com- 
ponent, monthly values for PV panels tilt angle, the height of wind 
turbine towers along with the PMS parameters ensure the capabil- 
ity of the proposed method to achieve the aim of optimization. In 
order to focus on the significance of PMS and show the influence of 
strategy variations on the performance and objectives of the sys- 
tem, a comparison is presented, too. In addition, determining the 
PMS parameters in the optimization procedure is confirmed as an 
efficacious concept. 

In the next section, the description of the system under study 
and the model of components are presented. Section 3 illustrates 
the problem statement including objective functions, constraints 
and the proposed method for the power management. In Section 
4, DEA and fuzzy multi-objective technique are explained. Finally, 
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Fig. 1. Configuration of stand-alone hybrid energy system under study. 


Fig. 2. Representation of solar elevation angle and PV panel tilt angle. 


the results and conclusions are represented in two subsequent 
sections. 


2. Model of system components 


In this work, a stand-alone hybrid energy system comprising 
wind turbines, PV panels, fuel cells, electrolyzers, hydrogen tanks, 
batteries and diesel generators is studied. A simplified diagram of 
the system is depicted in Fig. 1. The available power from renew- 
able sources is directly delivered to the load in order to provide the 
load demand. The excess or deficient power from RES is saved in or 
supplied from other equipments in the system based on the pro- 
posed dispatch strategy explained in Section 3.2. The model of each 
system component is briefly described in the following subsections. 


2.1. Wind turbine 


The output power of the wind turbine, Py, is calculated on the 
basis of its power curve given by the manufacturer. The effect of 
the wind turbine installation height and the roughness of the land 
surface of the towers on the output power is also considered [2,6]. 


2.2. PV panel 


To achieve the maximum output power from a PV panel, it is 
essential to consider the influence of the solar position and panel tilt 
angle. In [16], itis shown that such effects can change the total avail- 
able power up to 9.94% of the total maximum power. The panels 
are assumed to have rotation on one axis. The following equations 
lead to the incident irradiation on a tilted panel. Eq. (1) is employed 
to calculate the longitude at the equator (ô), with respect to solar 
ecliptic longitude and latitude (A, 0) and earth’s axis inclination to 
the plane of its orbit ([=23.439°). 


Bec - Pyre + Arc : Pfc 
Consrc = 


P : 
Brc - Pyre + Arc - Pfc (1 + Fef (— = Pray #t)) otherwise 


ô= sin™! (sin 6 cos I + cos @ sin A sin I) (1) 


The solar elevation angle (o), which is the angle between the direc- 
tion of the sun and the horizon is then estimated, as follows: 


o = sin”! [sin gy sin ô+ cos g cos 6 cos(LMST — À )] (2) 


where LMST stands for Local Mean Sidereal Time and g is the geog- 
raphy of the longitude [17]. The diffuse and reflected radiations are 
neglected. Using Fig. 2, the incident radiation on the tilted surface 
(Gi) can be expressed in terms of horizontal component of solar 
irradiation (G), as follows: 


Gn 
= sin o (3) 
Gp = Gsin(o + £) (4) 


i 


where Gp is the effective component of the solar radiation perpen- 
dicular to the panel. Based on the calculated Gp, the power available 
from PV panels at time step t considering the temperature effect is 
determined using the following equations: 


NCOT — 20 


Te(t) = Tat) + SSP Gace, 6) (5) 
Isc(t) = Use sre + Ki(Te(t) — 25°)] 20) (6) 
Voc(t) = Voc,stc — Kv - Tc(t) (7) 
P(t) = Ns Np Voc «Isc(t) - FF(E) (8) 


where NCOT is the nominal cell operating temperature (°C), T,(t) 
is the ambient temperature at t, T-(t) is the cell temperature at t, 
Vocstc and Iscstc are the module open-circuit voltage and short- 
circuit current under STC, with K; and Ky as their corresponding 
temperature coefficients; Pyy(t) is the power of array with Np mod- 
ules in parallel and N; modules in series and FF(t) is the fill factor 
[6]. 


2.3. Fuel cell 


The formulations, which represent the balance between the 
hydrogen as input and the electric power as output in steady state 
condition, is the required fuel cell model used in the design pro- 
cedure. However, for the sake of ensuring maximal efficiency of 
the system and preventing degradation of fuel cells due to abun- 
dance of start-up and shut-down cycles, a hysteresis band is also 
embedded into the model, as in [12]. 

The fuel cell hydrogen consumption, Consrc (kg/h), relative to 
its output electric power is computed by the following equation 
[15]: 


for P/PN Fc < Pmax ef 


(9) 


where Prc is the output power of fuel cell in kW, Py rc is the nomi- 
nal output power, Arc and Bre (kg/kWh) are the coefficients of the 


1580 S. Abedi et al. / Renewable and Sustainable Energy Reviews 16 (2012) 1577-1587 


consumption curve. Pmax ef, in terms of percentage of Py rc, is the 
output power, at which the efficiency of the fuel cell is maxi- 
mum and Fe is a factor to consider the high consumption above 
Pmax ef- The hydrogen consumption parameters for all fuel cells are 
assumed to be Agc=0.05 and Brc = 0.004 kg/kWh, Pmax ef =0.2 and 
Fop=1 [15]. By using these values, the fuel cell efficiency is 46% at 
the maximum and 31% at rated powers. 


2.4. Electrolyzer and hydrogen tank 


The input electrical energy dependence on the hydrogen mass 
flow is modeled, as follows [15]: 


Conse = Be - Que + AE: Q (10) 


where Qy: is the nominal hydrogen mass flow (kg/h), Q is the 
hydrogen mass flow (kg/h), Ag and Bg are the coefficients of the 
consumption curve in kW/kg/h. The parameters of the electrical 
energy consumption for all electrolyzers are assumed to be the 
same as parameters given in [15], i.e., Ag=40 and Be =20 both in 
kW/kg/h. 

The maximum hydrogen tank capacity is assumed to be equal to 
10kg and the total capacity of the hydrogen storage is determined 
by the design program as the number of hydrogen tanks (Nang). 
The minimum allowed level of the hydrogen (in tanks) is also a 
parameter determined by the proposed PMS. 


2.5. Battery 


The battery bank with the total nominal capacity of Cn (Ah), can 
serve energy to the load until the maximum depth of discharge 
(DOD) or SOCpin is reached. SOCmin Of battery bank along with 
SOCmax are taken into account as two of the parameters of PMS, 
which should be optimized by the developed program and should 
not exceed the values introduced by the manufacturer. 

The SOC of battery bank, in each simulation time step, is calcu- 
lated by applying the following equation: 


SOC(t) = SOC(t — 1) + ng 22 . 100 (11) 
where SOC(t) is the batteries’ SOC at time step t, ng is the battery 
round-trip efficiency and Pg(t) is the power charged in or dis- 
charged from battery bank at time step t [6]. ng is approximately 
equal to 80% in charging and 100% in discharging modes [6]. Pg(t) 
is determined due to the mismatch of power between load and 
RES and the share associated with batteries to supply the defi- 
cient power or to save the excess power. The mentioned share is a 
parameter of PMS and optimized in the proposed algorithm. 

The number of series connected batteries depends on the DC bus 
voltage (Vgys) and the nominal voltage of each battery, i.e.: 


— Vsus (12) 


The total number of batteries in battery bank is equal to ne multi- 
plied by nb (the number of batteries in parallel), which should be 
determined in the optimization algorithm. 

In order to model the effect of the number of charge/discharge 
cycles on the lifetime of batteries, the model used in HOMER soft- 
ware [22] is used. This model is based on the battery “cycles to 
failure” curve, which is more explained in detail in [11]. 


2.6. Diesel generator 


The fuel consumption of the diesel generator, Consg, in terms of 
its output power is written, as follows: 


Consg = Bg - Py c + Ag: Pc (13) 


where Ag and Bg (1/kW) are fuel consumption curve coefficients 
provided by the manufacturer, Pg (kW) is the output power and 
Pyg (KW) is the nominal output power of the diesel generator. 

The values assigned to Aç and Bg are 0.246 and 0.08145 1/kWh, 
respectively, for all diesel generators, which have been used in [15], 
too. 


3. Problem statement 
3.1. Objective functions and constraints 


The objective functions, considered in the optimization prob- 
lem, are as follows: 


e The overall cost of the system discounted to the year of 
installation (NPC) including the investment cost, operation and 
maintenance cost during the lifespan of the system, replacement 
cost, and fuel cost of diesel generators. All of the cost terms asso- 
ciated with the system components are based on the data in [15]. 

e The total fuel emissions produced by diesel generators during the 
total lifespan of the system. Thanks to the employed fuzzy multi- 
objective optimization technique, the model implemented for the 
fuel emission is adequately accurate and there is no requirement 
for taking it into account in terms of cost. As the most significant 
gas existing in the diesel generator exhaust is CO2, which also 
causes green house effects, the produced CO% is considered to 
represent the fuel emission. 

e The total Loss of Load Probability (LPSP), which represents the 
total unmet energy. This index is the most commonly used index 
in related works [2]. LPSP is proportional to the unmet load and 
is written, as follows: 


psp = LOE (14) 


n-1D(h) 


where D(h) is the load demand in kWh at time step h, and LOEE, 
standing for Loss of Energy Expected, is defined by the following 
equation: 


LOEE = YO FILOE(h) (15) 
E[LOE] = $ “Q(s)- f(s) (16) 
seS 


where Q(s) represents the amount of the unmet energy when the 
system experiences the state s and f(s) is the probability associated 
with the occurrence of the state s. 

The control variables, which should be used in simulations to 
calculate the objective functions, as well as their boundaries, are 
listed, as follows: 


- A k x 2 matrix [[Nx], [T;,.]] of size parameters of the system com- 
ponents shown in Fig. 1, where, N; is the number and T;, is the 
type of component k. Both are positive integers. 

- Avector of installation parameters including monthly tilt angles 
of PV panels (0° < 6, < 90°, k=1,...,12) and wind turbine tower 
height (h > 0), 

- A vector of PMS parameters (described in the coming section) 
including: 

- Monthly charge shares associated with n—1 types of storage 
media (0<CS; <1, where, n is the total types of the storage 
media existing in the system), 

- Monthly discharge shares associated with m — 1 types of backup 
devices (0 < DS; <1, where, m is the total number of backup 
devices capable of feeding power to the load), 
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- Upper and lower limits on the stored energy level 
of all the storage media (SOCjminrated <SOGmin < 0-5, 
0.5 < SOC, max < SOC;,max,rated, 0 < SOC < 1). 


In addition to the boundaries that control variables should not 
take values beyond them, there are operational constraints to be 
mentioned as follows: 


- The specified maximum and minimum SOC of batteries (SOCpin, 
SOCmax), 

- The rated current or power of all devices, 

- The minimum power injected to the electrolyzers and the min- 
imum supplied power from fuel cells and diesel generators, at 
which they can start-up. This constraint is considered to avoid 
the multitude of the start-up/shut-down cycles or short periods 
of frequent utilization of these devices, and enhance their life- 
time. This operation limit has been introduced as Hysteresis Band 
Range (HBR) in [12]. Other imposed constraints are illustrated in 
the next section. 


3.2. Proposed method for power management of system 


In general, the PMS is so important that the operation, relia- 
bility, cost and lifetime of the system is affected with even minor 
alterations in the strategy. In the proposed dispatch strategy, to effi- 
ciently utilize the RES, the available power from RES, i.e., Pres(t), is 
directly delivered to the load in order to provide the load demand, 
Pigaq(t). Other equipments are considered without any predefined 
priorities and there is not any preference to some components. The 
priorities are determined by the optimization process. The excess 
or deficient power of RES, i.e., Prem(t), is saved in or supplied from 
other equipments in the system, i.e.: 


Prem(t) = Pres(t)- Pioaa(t) (17) 


where Prem(t) might be a negative or positive value. 

For Prem(t) > 0: In this case, Pres(t) is more than the load demand 
at time step t, and the remaining power should be delivered to the 
storage devices namely electrolyzers and battery bank. Therefore, 
based on the proposed method, the optimization procedure should 
assign Charge Shares (CS) for battery bank and electrolyzers so that 
Prem(t) is optimally shared among them considering constraints. 
The power, to be stored in storage device J, i.e., P;(t), is written, as 
follows: 


P;(t) = CSi(k) - Prem(t) i=1,...,n (18) 


where n is the total number of storage media and CS;(k) is the CS 
of storage device i in month k (k=1,...,12). It should be mentioned 
that all CSs fall within the range (0,1). 

Now, the status of all storage devices is checked one after 
another (as seen in Fig. 3), considering their constraints and lim- 
itations on the storable power. If any constraint is violated, P;(t) 
of the device should be modified so that the violated constraint is 
satisfied. These constraints violations or incompatibilities and the 
necessary modification to P;(t) in each case are, as follows: 

Case (1-1): The power assigned to device i, P;(t), exceeds its rated 
power. In this case, it is set to the rated power (P,atea,i) of the device, 
i.e.: 

P¥(t) = Prated,i (19) 


Case (1-2): P;(t) exceeds the remaining storage capacity of device i. 
In this case, the power is injected to the device until the maximum 
allowable level of the stored energy is reached (SOCmax): 


Prated,i 


Ps(t) = 


(SOCmax,i — SOC(t — 1)) (20) 


The value of SOCmax is obtained within the context of the optimiza- 
tion procedure. 

Case (1-3): P;(t)is beyond the HBR. Hence, the power to be drawn 

from the device does not fall within the allowable range and no 
power can be delivered, i.e.: 
P*(t)=0 (21) 
Case (1-4): The operating limitations of the device considering its 
dynamics are not satisfied. For instance, the change in the power 
share of the device from time step t to {+1 is more than the maxi- 
mum tolerable rate of the power change due to the system inertia 
(i.e., the power ramp rate, represented by r). This inertia implies the 
impossibility of timely tracking very rapid input power fluctuations 
[23]. 


P*(t)=r- At (22) 


In above mentioned cases, the power assigned to other devices, 
P;(t), should also be modified because the constraint did not allow 
the device i to fulfill the storage of its initially assigned P;(t), calcu- 
lated in (18). Therefore, a portion of P;(t) is still remained, namely 
Premi(t), which other devices may be able to absorb it. The modifi- 
cation procedure is named Charge Share Correction (CSC), written 
as follows: 


Prem,i(t) = Pi(t) — P#(t) (23) 
CS*(t) = 0 (24) 


CS; 
1+ a -5 
et kei k 
ii 


P*(t) = Pi(t) + C5*(t)- S Prem xlt) (26) 
k=1 


so-s ( forj#i,j=1,...,n (25) 


where CS; is the charge share of all devices except device i, and P(t) 
is the modified power assumed for the storage device j. Eq. (24) 
means that the device i has absorbed all possible amount of power 
it could, and there is no more power stored in it during the current 
time step t. If one of the constraints corresponding to this device is 
violated, the same modification procedure is then performed, for 
that constraint. 

For Pyem(t) < 0: In this case, the load demand is partially supplied 
by RES. Hence, Discharge Shares (DS) should be determined and 
assigned to the back-up devices. Utilizing these back-up devices, 
the effort is made to fully compensate the remaining power by the 
power from suppliers other than the RES (Ppackup ), and try to achieve 
Phackup = |Prem|. 

The power to be drawn from backup device i, P;(t), is written as 
follows: 

P,(t) = DS;(k) - Prem(t) i=1,...,m (27) 
where m is the total number of back-up devices capable of feeding 
power to the load and DS;(k) is the discharge share of the back-up 
device iin month k (k=1,...,12). All DSs fall within the range (0,1). 

As shown in Fig. 3, the status of all back-up devices is checked 
one after another considering their constraints and limitations on 
the power that they are able to serve in the current time step. If 
any constraint is violated, P;(t) of the device should be modified, as 
follows, so that the violated constraint is satisfied. 

Case (2-1): P(t) exceeds the rated power of device i. (Eq. (21) 
should be applied.) 

Case (2-2): The stored energy in the device in the current sim- 
ulation time step is not as adequate as to supply P,(t). In this case, 
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Case (1-1) Yes 
occurred 
Eq. (19), CSC 


Case (1-2) 
occurred 


Case (1-3) 
occurred 


Y Case (1-4) 


es 
occurred 
CSC 


Yes Case (2-1) 
occurred 
Eq. (19), DSC 
No 
Case (2.2) Yes 


occurred 


No 
Yes Case (2-3) 
occurred 
Eq. (21), DSC 


No 


Yes 
Eq. (22), DSC 
Calculate supplied and 
unsupplied energy 
<a 
Yes 


Case (2-4) 
occurred 


No 


Fig. 3. Flowchart of proposed PMS. 


the power is drawn from the device until the minimum allowable 
level of the stored energy is reached (SOCmin,i): 


P f 
T (SOC;(t — 1) 
The value of SOCmin is obtained within the context of the opti- 
mization procedure. 
Cases (2-3) and (2-4): These cases are the same as cases (1-3) 
and (1-4), respectively. 
In the above mentioned cases, the power assigned to other 
devices, P;(t), should also be modified because the constraint did 


P*(t) = SOCmnin,i) (28) 


not allow device i to fulfill the supply of its initially assigned P;,(t), 
calculated in (27), and a portion of P;(t) is still remained, namely 
Premi(t), which other devices may be able to supply it. The modifica- 
tion procedure is named Discharge Share Correction (DSC), written 
the same as CSC (Eqs. (23)-(25)), with only replacing all CS variables 
by DS, and n by m. 

While any of the above situations take place and CSC or DSC is 
performed, it is possible that again, the modified dispatch shares are 
incompatible with system operation constraints. In this case, the 
CSC or DSC is performed for the second time so that the constraints 
are satisfied. This is the last effort for the determination of dispatch 
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Fig. 4. Membership functions of objective functions and constraints. 


shares. Then the remaining excess or deficient power is considered 
as discarded or unmet power in the current simulation time step, 
respectively (as presented in Fig. 3). 


4. Optimization algorithm 


The multi-objective optimization algorithm proposed in this 
paper, takes the advantage of the differential evolution algorithm 
accompanied with fuzzy-multi objective technique. 


4.1. Multi-objective optimization using fuzzy technique 


The multi-objective problem is generally solved by three types 
of methods. The first one is pareto-based approach to get a set 
of non-dominated solutions in the process of optimization [24]. 
The second one is the coefficient method and the last method is 
transforming the multiple objective function into a single objec- 
tive model and optimizing it through single objective strategies 
[24]. In this method, the determination of appropriate values for 
the coefficients is of very important. In addition, the results are 
highly dependent on and sensitive to selected values for coeffi- 
cients. In the method developed by Bellman and Zadeh [25], the 
single objective problem is achieved by maximizing the minimum 
degree of satisfaction among membership functions. 

The basic idea in fuzzy optimization is to optimize objective 
function and constraints simultaneously [26]. The fuzzy decision 
is marked out due to the intersection of fuzzy objectives and fuzzy 
constraints. The first operation is the fuzzification process of the 
merged objective function and the constraints. In this procedure, 
two types of function u(x) are defined for each objective function 
and constraint, as shown in Fig. 4. In this figure, the minimum value 
for each objective is obtained by the single objective optimization 
and the maximum value is specified by the initial set of control 
variables in the optimization algorithm. According to this method, 
it is possible to change the effectiveness of any objective function 
by reducing or increasing its specified maximum value. In other 
words, when the specified maximum value of an objective function 
decreases, it is considered to be more important in the optimiza- 
tion than the previous state and vice versa. These membership 


functions are initially combined by “and” operator (minimum). This 
procedure can be expressed by the following equation: 


Mp(x) = min(ufı(x), uf2(x), .. -, Mer (X), Mca(X), ++) (29) 


where up(x) represents the membership function of the optimal 
decision function and is treated as the evaluation value in the opti- 
mization algorithm. 

The membership values express the degree of satisfaction for 
each objective. Highly satisfied objective is given a low value, 
though lowly satisfied one is assigned a high value. Hence, the 
multi-objective problem can be transformed into the following 
maximization problem subject to a crisp constraint set: 


Max up(x, u)s.t. H(x,u)=0, C(x,u)<0 (30) 


where up(x,u) is called the fuzzy index. 
4.2. Differential evolution algorithm 


Differential Evolution Algorithm (DEA), introduced by Storn 
and Price [27], is a simple population based stochastic parallel 
search evolution algorithm for global optimization and is capable 
of handling non-differentiable, nonlinear and multi-modal objec- 
tive functions [28]. In DEA, the population consists of real-valued 
vectors with dimension D, which is equal to the number of control 
variables. The initial population with the size Np, is uniformly dis- 
tributed in the search space, falling within variables’ boundaries. 
The procedure of this algorithm is shown in Fig. 5. The algorithm is 
described in the following steps: 

Step (1), Initialization 

For each control variable k with lower bound ae and upper 


bound x% ax initial values are randomly and uniformly chosen in 


i k ky. 
the interval [x0 in Xmax]: 


k k 


Fo = XÉ in + rand [0,1] x (xhax —%4;,) iell, Np], 


x mii ke[1, D](31) 


Step (2), Evaluation 

In DEA, after generation of the initial population Xg, all of its 
consisting vectors are evaluated by the calculation of the objective 
function in its subprogram, and then the vectors are sorted accord- 
ing to their objective function value. In the present work, DEA is 
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Fig. 5. Procedure of the used DEA. 


accompanied with fuzzy technique to handle the multi-objective 
problem. Hence, this step has some differences, described in Section 
4.4. 

Step (3), Iteration 

This step is the main procedure of DFA to conduct the final 
results, which is performed iteratively, until either an acceptable 
solution has been reached, or for at least a specified number of iter- 
ations, no further improvement in the solution is obtained, or a 
predefined number of iterations are completed [28]. 

Step (3-1), Mutation: 

In DEA, the perturbations applied to the current population to 
generate new population are derived from the current population 
itself, and no predefined probability density function is considered. 
For a given vector of control variables, X; ç, in the population, two 
vectors (X,, ¢ and X,, ¢) are randomly selected such that the indices 
i, r2 and r3 are mutually different. Then, the weighted difference 
between two vectors is added to the vector with the best fitness 
function, Xp g, written as follows: 


Vict = Xp,G +F- (Xn,G — Xr,G) 12,73 €{1,2,..., Np} (32) 


where F, a constant from (0, 2), characterizes the amount of the 
movement of vectors within the search space. 

Step (3-2), Crossover: 

In this step, the trial vector U; +1 is developed from the elements 
of the target vector X; ç and the elements of the mutant vector V;¢+1. 


The elements of the mutant vector enter the trial vector with prob- 
ability CR as follows: 


ee { ViiG+1 if rand; ; < CR orj = land (33) 

na Xj if randj; > CRandj + Irana 
where randjie [0,1], Irana is an integer randomly selected from 
{1,2,..., D]. Iranq ensures that one of the trial vectors is selected 
among the mutant vectors. 

Step (3-3), Modification and constraints verification: 

This step prepares the trial vectors for the next step by veri- 
fying the variables whether they meet their constraints, given in 
Sections 3.2 and 4.4. Otherwise, for the vectors containing vari- 
ables with unsatisfied constraints, mutation and crossover steps 
are repeated for a specified number of times or until the problem 
is fixed. Nevertheless, these vectors will be regenerated, similar to 
step (1). 

Step (3-4), Selection: 

The vectors generated in the previous step are evaluated and 
sorted as described in the step (2). Then, the trial vector Ujc¢+1 is 
compared with the corresponding vector X;¢ and the one with the 
better fitness value is admitted to the next iteration as Xj,¢+1: 


_ J Uic+ iff(Uic+1) <f(Xi,e) ; 
Xic = E otherwise ie[1, Np] (34) 


4.3. Representation of uncertainty 


Natural characteristics of wind and solar energy impose uncer- 
tainty in their design and operation. Hence, considering various 
possible scenarios in the model of these resources can lead to 
more realistic decisions. The uncertain parameters are expressed 
by probability distributions, showing the range of values that a 
parameter could take, and also accounting for the probability of the 
occurrence of each value in the considered range [12]. In this paper, 
Weibull and Beta PDFs are used for the wind speed and the solar 
irradiation, respectively [21]. Using this solution, seven scenarios 
are generated, to model the uncertainty of these resources. 


4.4. Implementation 


A program in MATLAB environment is developed to determine 
the best set of design and PMS parameters. Furthermore, the con- 
vergence and computational time for the algorithm is evaluated. 
The same system is simulated via two other optimization methods 
widely applied to renewable energy [29], namely Genetic Algo- 
rithm (GA) [30] and Linearly Decreasing Inertia Particle Swarm 
Optimization (LDI-PSO) [31,32] and their corresponding results are 
compared in the Section 5.2. The performance of the proposed PMS 
is also evaluated by comparison with another PMS. 

There are two main modules in the developed program, the 
optimization module and the system simulation module. The first 
module generates and modifies a valid vector of control variables 
to be evaluated by the system simulation module in each iteration. 
When the latter module is invoked, the input vector is processed. 
The module simulates the operation of the system for each scenario, 
calculates the operating point of each system component and deter- 
mines /4¢;(x) and uc(x), and the resultant up(x), using the resultant 
values for objective functions. These values are transformed into a 
single value using the described fuzzy technique and become the 
evaluation value for the considered scenario. The final evaluation 
value for an input vector is the average of values obtained by using 
data of all generated scenarios [12], which is then admitted back to 
the optimization module. 

Therefore, a valid control vector is obtained by assigning val- 
ues to the control variables that lie within allowed bounds of each 
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Fig. 6. LOEE for (a) PMS (a) and (b) PMS (b). 


individual. However, some modifications to the original optimiza- 
tion algorithm are made to match the specific characteristics of the 
present problem, as follows: 


- Some of the control variables, including number and type of each 
component, must be non-negative integers. Therefore, the opti- 
mization procedure should handle the discrete variables in all 
initialization, crossover and mutation steps. 

- The upper limit, associated with CS; and DS; (i=2,...,n) float in 
accordance with the values of CS; and DS; (k=1,...,i— 1), because 
CS; and DS; must satisfy the equality constraint are = 


ye, = 1. This constraint modification is applied to the value 
that the optimization procedure assign to CS; and DS;, from device 
i=1 through n, one by one, i.e.: 


i-1 
csmax 1, Cpr} "cs, S39. Gh) (35) 
k=1 


A similar equation is used by replacing CS by DS and n by m. 

The integration of the fuzzy model with the optimization pro- 
cedure for handling the multi-objective problem consists of 
developing membership functions of objectives as well as con- 
straints imposed on the operation of the system. By solving 
a single objective optimization for each objective function, its 
minimum and maximum values are obtained and the corre- 
sponding membership function is achieved. By the way, using 
maximum and minimum power ramp rates of each compo- 
nent (ramp rate; pin and ramp rate; max), rated limitations on SOC; 
(SOC; min and SOC; max) and power boundaries of each component 
(rated poweri min and rated power; max), all as the operation con- 
straints, their corresponding membership functions are obtained 
(as shown in Fig. 4). 

To treat the uncertain parameters, several scenarios are generated 
and simulated. Hence, all the three objective functions introduced 
in Section 3.1 are calculated for all the scenarios. The algorithm 
utilizes the concept of the “expected objective function”, which 
aggregates all the objective function values corresponding to 
the simulated scenarios representing the uncertain parameters, 
in the form of an average value [12]. Then, using the expected 
objective function values, the fuzzy index is calculated, and the 
algorithm continues with the same algorithmic operations found 
in DEA. In this way, the results represent the solution with an 
overall consideration of all scenarios associated with the uncer- 
tain parameters in the system. 


5. Case study and results 
5.1. System data 


The hybrid energy system under study has been described in 
Section 2. In practical system planning, the designer has to choose 
the system components from a set of commercially available ones. 
As acase study, various models and types of components with dif- 
ferent rated powers, costs and other specification [15], are the input 
to the design procedure. The system location is a region in Arde- 
bil city (38°14’57”N/48°18/5’’E), in the north west of Iran. For this 
region, the hourly data for solar irradiance, wind speed and ambient 
temperature all provided within one year period [2]. The consid- 
ered load profile is the IEEE reliability test system [33] with 30 kw 
peak power. 


5.2. Results 


The data are analyzed to determine the optimum design and 
PMS parameters for the case study. The results of the used opti- 
mization algorithms, i.e., GA, LDI-PSO and DEA, are presented in 
Table 1. It can be seen in this table that DEA revealed the best per- 
formance with respect to optimization of objective functions and 
the resultant fuzzy index. Hence, DEA is nominated to demonstrate 
the results of this problem in more detail. Since, this problem should 
be solved in off-line mode, the computation time and convergence 
speed are not of great concern. 

In order to well present the applicability and efficacy of the pro- 
posed PMS, which is the significant idea of this paper, a simulation 
is also carried out by using another PMS, referred to as PMS (b). 
Based on the PMS (b) which is analogous with the PMS (dispatch 
strategy) employed in the HOMER software [8], Prem(t) in time step 
t is dispatched among components by the priorities defined by 
designer, unless the marginal cost of energy storage or production 
of the dispatched power by the components is in contrast with 


Table 1 
Results of three evolutionary algorithms. 
PMS (a) PMS (b) 
GA LDI-PSO DEA DEA 
Cost (€) 1,042,242 967,375 912,560 918,337 
Fuel emissions (kg) 2086.2 2632.8 1827.7 2888.6 
LPSP (W/h) 0.0067 0.0083 0.0072 0.0075 
Number of iterations to converge 29550 196 227 153 
Computation time (s) 173,520 121,311 134548 72,549 
Fuzzy index 0.9980 0.9973 0.9996 0.9915 
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Table 2 
Optimized values for design parameters of the case study. 
Design data DEA 
Equipments data 
Number 9 
PV panel Type 10 
g i Number 6 
Wind turbine Type 6 
Number 28 
Battery Type 5 
Number 6 
Fuel cell Type 4 
Number 8 
Electrolyzer Type 4 
Diesel generator ah da i 
Hydrogen tank Number 6 
Inverter nominal power (kVA) Type 35 
Installation data 
Wind tower height (m) 32 
67 
58 
63 
46 
38 
PV panels monthly 31 
tilt angle By — B12 (°) 22 
30 
26 
43 
59 
53 


the predefined priorities. In this case, the priority of components 
to store or produce Prem(t) is sorted ascendingly according to 
their corresponding cost of energy storage or production. The 
optimization results using PMS (b) are included in Table 1, too. It 
can be observed that PMS (a) have better results than PMS (b). The 
fuel emission is comparably less than that of PMS (b), showing 
that the utilization of back-up and storage devices other than the 
diesel generator is more effective when PMS (a) is employed. Fig. 6 
presents the LOEE in both PMS cases, confirming the better LPSP of 
Table 1 for the case of PMS (a). 

Table 2 presents the result of design parameters, including the 
type and number of each component. The installation data con- 
sisting of the monthly PV panels tilt angle and wind turbines 
installation height are also listed in the table. 

The monthly values of charge and discharge shares are pre- 
sented in Figs. 7 and 8, respectively. Fig. 7 presents that the DS 
of fuel cells is more significant in winter’s months than summer’s 
months and the opposite is true for diesel generators. The values 
obtained for DS parameters also have great influence on the CS 


—+t— Battery 
4 = Diesel 


—*— Fuel cell 


Fig. 7. Monthly dispatch shares of backup devices for providing the deficient 
demand. 
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—— Battery 


—*- Electrolyzer 


Fig. 8. Monthly dispatch shares storage devices for saving the excess power. 


Table 3 

Optimized values of limitations of storage devices. 
Storage device limitations DEA 
Minimum state of charge of batteries (SOCmin %) 21 
Maximum state of charge of batteries (SOCmax %) 100 
Minimum level of Hydrogen in tanks (%) 33 


values. In fact, CS and DS parameters are interdependent for compo- 
nents that utilize the same energy storage device in both charging 
and discharging states. For instance, electrolyzers and fuel cells use 
the same energy container namely hydrogen tanks to store energy 
and consume it, respectively. Table 3 contains optimized values of 
operating limitations of storage devices. 

The attained values for parameters of PMS by represent nearly 
the overall optimum operating condition of the system considering 
a compromise among all the considered objectives. The parameters 
are separately calculated for each month, representing the influ- 
ence of monthly and seasonal changes in load demand and climatic 
patterns on the utilization of the hybrid system. Hence, the results 
are adapted with any circumstances that alter the operation of the 
system resulting in more accuracy and optimality. 


6. Conclusions 


This study has been dedicated to review optimal design and 
operation strategy selection of hybrid RES-based stand-alone 
energy systems, including various generators and storage media. 
Considering resource uncertainties, associated with wind speed 
and solar irradiation, a novel method has been proposed, to deter- 
mine the power management strategy of the system along with 
sizing parameters of system components so that the overall cost 
of the system during its useful lifetime, unmet load and pollutant 
emissions have been minimized. The PMS parameters are monthly 
charge shares (CS) of storage devices and monthly discharge shares 
(DS) of generator devices, limitations on their start-up/shut-down 
power thresholds (HBR) and on their level of available energy (SOC). 

Based on the proposed method, the values of PMS parameters, 
the sizing and design parameters have been determined by the 
iterative optimization algorithm. The PMS parameters have been 
separately determined for each month to adapt monthly variations 
in the load and climate patterns. Then in each iteration, these values 
have been applied to the hourly simulation of the system opera- 
tion and evaluated in each single time step to meet the operational 
constraints of the components, consisting of the nominal power, 
SOC min and SOCmax, HBR, and power ramp rate of each device. Oth- 
erwise, these values have been modified by either of two operators, 
namely CSC or DSC, in that time step. These modifications have 
ensured that while the optimum operation strategy and commit- 
ment of components has been attained, all constraints associated 
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with operational characteristics of the components have been sat- 
isfied. However, if after completion of all steps of simulation in a 
time step, the constraints were still unsatisfied, some discarded or 
unsupplied energy remains. The summation of these values during 
all time steps of the simulation has been involved in the objective 
functions to be minimized. 

The proposed method has tested on a hybrid 
PV-wind-diesel-hydrogen-battery system. To efficiently handle 
the nonlinear mixed-integer multi-objective optimization prob- 
lem, a modified version of DEA accompanied by fuzzy technique 
has implemented. The system has been optimized for summations 
of all objective values of each function calculated for all scenarios 
generated by applied uncertainty models. The presented results of 
the solution can be categorized into three groups: 


- Size (design) data: The optimum number and type of each com- 
ponent of the system to be installed. 

- Installation data: The optimum values for monthly inclination 
angles of solar panels and the optimum tower height for wind 
turbines installation to maximize the exploitable power from RES. 

- Operation strategy (PMS): The monthly CS and DS, and boundaries 
associated with SOC. 


By comparing the results of design parameters and objective 
functions values (accompanied by that of the system operation 
and unsupplied energy), with the results of applying a user- 
defined PMS, it has been concluded that the design decisions and 
operation of such systems are highly affected by the employed 
PMS; and the superiority of the proposed PMS has been con- 
firmed. It should be mentioned that the proposed method can 
be effectively employed for any composition of hybrid energy 
systems. 
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